home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Shareware Grab Bag
/
Shareware Grab Bag.iso
/
090
/
byt1186b.arc
/
APPROX.LBR
/
REMES.C
< prev
next >
Wrap
Text File
|
1986-04-11
|
6KB
|
272 lines
/* remes.c
*
* This is an interactive program that computes least maximum polynomial
* and rational approximations.
*/
#define P 15 /* max total degree of polynomials, + 2 */
#define N 20 /* number of items to tabulate for display */
extern double PI; /* 3.14159... */
static int IPS[P] = {0,}; /* simq() work vector */
static double AA[P*P] = {0.0,}; /* coefficient matrix */
static double BB[P] = {0.0,}; /* right hand side vector */
static double param[P] = {0.0,}; /* solution vector */
static double xx[P] = {0.0,}; /* points in approximation interval*/
static double ref[N+1] = {0.0,}; /* function values for display */
static int n = 0; /* degree of numerator polynomial*/
static int d = 0; /* degree of denominator polynomial*/
static int nd1 = 0; /* n + d + 1 */
main()
{
int i, ip, j, ncheb, neq, relerr;
double a, apstrt, apwidt, b, c, x, y, z;
char s[40];
double abs(), approx(), cos(), func();
int simq();
printf( "\nRational Approximation by Remes Algorithm\n\n" );
START:
printf( "Relative error (y or n) ? " ); /* Ask for error criterion */
gets( s ); /* read in a line of characters */
relerr = 0;
if( s[0] == 'y' )
relerr = 1;
printf( "Degree of numerator polynomial? " );
gets( s ); /* read line */
sscanf( s, "%d", &n ); /* decode characters */
printf( "Degree of denominator polynomial? " );
gets( s );
sscanf( s, "%d", &d );
printf ( "Start of approximation interval ? " );
gets( s );
sscanf( s, "%E", &apstrt );
printf ( "Width of approximation interval ? " );
gets( s );
sscanf( s, "%E", &apwidt );
nd1 = n + d + 1;
/* remes.c 2 */
/* Supply initial guesses for points in approximation interval */
if( d == 0 )
{ /* there is no denominator polynomial */
neq = n + 2; /* The number of equations to solve */
ncheb = n + 1; /* Degree of Chebyshev error estimate */
/* Extrema of Chebyshev polynomial */
a = ncheb;
for( i=0; i<neq; i++ )
{
xx[i] = apstrt
+ 0.5 * apwidt * (1.0 - cos( (PI * i) / a ) );
}
}
else
{ /* there is a denominator polynomial */
neq = nd1;
ncheb = nd1;
/* Zeros of Chebyshev polynomial */
a = 2.0 * ncheb;
for( i=0; i<ncheb; i++ )
{
xx[i] = apstrt
+ 0.5 * apwidt * (1.0 - cos( PI * (2*i+1) / a ) );
}
c = 0.0; /* deviation at solution points */
}
/* calculate function table for error curve display */
a = apwidt/N;
b = apstrt;
for( i=0; i<=N; i++ )
{
ref[i] = func(b); /* func is the f(x) to be approximated */
b += a;
}
/* remes.c 3 */
LOOP:
/* Display old values of guesses and let user change them if desired */
/* First do the deviation guess if rational form */
if( d > 0 )
{ /* there is a denominator */
c = abs(c); /* deviation at solution points */
printf( "deviation = %.4E ? ", c );
gets( s );
if( s[0] != '\0' ) /* if input is not a null line, */
sscanf( s, "%E", &c ); /* then decode the number */
}
else /* no denominator: */
c = 1.0; /* coefficient of the deviation */
/* Read in guesses for locations of solution */
for( i=0; i<neq; i++ )
{
printf( "x[%d] = %.4E ? ", i, xx[i] );
gets( s );
if( s[0] != '\0' )
{
sscanf( s, "%E", &x );
xx[i] = x;
}
}
/* remes.c 4 */
/* Set up the equations for solution by simq() */
for( i=0; i<neq; i++ )
{
ip = neq * i; /* offset to 1st element of this row of matrix*/
x = xx[i]; /* the guess for this row */
/* right-hand-side vector */
y = func(x); /* accurate function value f(x) */
if( d > 0 )
{ /* add the deviation if rational form */
if( relerr ) /* relative error criterion */
y = y * (1.0+c);
else /* absolute error criterion */
y = y + c;
}
/* insert powers of x[i] for numerator coefficients */
z = 1.0;
for( j=0; j<=n; j++ )
{
AA[ip+j] = z;
z = z * x;
}
/* insert denominator terms, if any */
if( d > 0 )
{
z = 1.0;
for( j=0; j<d; j++ )
{
AA[ip+n+1+j] = -y * z;
z = z * x;
}
BB[i] = y * z; /* right hand side vector */
}
else
{ /* no denominator */
BB[i] = y; /* right hand side vector */
y = c;
if( relerr )
y = y * BB[i];
AA[ip+n+1] = y;
}
c = -1.0 * c; /* switch sign of deviation for next row */
}
/* Solve the simultaneous linear equations */
simq( &AA[0], &BB[0], ¶m[0], neq, 0, &IPS[0] );
/* remes.c 5 */
/* Display the results */
j = 0; /* printout variable */
ip = 0; /* solution vector counter */
printf( "Numerator coefficients:\n" );
for( i=0; i<=n; i++, j++, ip++ )
{
if( j >= 7 )
{
printf( "\n" );
j = 0;
}
printf( "%.4E ", param[ip] );
}
if( d > 0 )
{
j = 0;
printf( "\nDenominator coefficients:\n" );
for( i=0; i<d; i++, j++, ip++ )
{
if( j >= 7 )
{
printf( "\n" );
j = 0;
}
printf( "%.4E ", param[ip] );
}
}
else
printf( "\nDeviation: %.4E", param[ip] );
/* Display table of function and approximation error */
printf( "\n\n x func approx error\n" );
a = apwidt/N;
b = apstrt;
for( i=0; i<=N; i++ )
{
x = b + i * a;
z = approx(x);
y = z - ref[i];
if( relerr && (z != 0.0) )
y = y/z;
printf( "%.3E %.3E %.3E %.3E\n", x, ref[i], z, y );
}
printf( "Another iteration (y or n)? " ); /* Ask what to do next */
gets( s );
if( s[0] == 'y' )
goto LOOP;
if( s[0] == 'x' )
exit(0);
else
goto START;
}
/* remes.c 6 */
/* This subroutine computes the rational form P(x)/Q(x)
* using coefficients from the solution vector param[]
*/
double approx(x)
double x;
{
double yn, yd;
int i;
yn = param[n]; /* highest order numerator coefficient */
for( i=n-1; i>=0; i-- ) /* work backwards toward the constant term */
yn = x * yn + param[i];
if( d > 0 )
{
yd = x + param[n+d]; /* highest order coefficient = 1.0 */
for( i=n+d-1; i>n; i-- )
yd = x * yd + param[i];
}
else
yd = 1.0; /* if there is no denominator */
return( yn/yd );
}
/* Put here an accurate routine for the function
* to be approximated
*/
double func(x)
double x;
{
double cos(), sqrt();
double y;
y = sqrt(x);
return( cos(PI * y / 2.0) );
}